\documentstyle[epic]{article}
\title{C++ Matrix Class}
\author{E. Robert Tisdale}
\begin{document}
\bibliographystyle{plain}
\maketitle

\begin{abstract}
This paper presents a simple, fast, efficient C++ Matrix class designed
for scientists and engineers.  The basic operations were implemented
in C++ using the Gnu \verb"g++" compiler (version~2.3.3) on Sun~4 computers
running UNIX (SunOS Release 4.1.1).  More sophisticated algorithms
were implemented by interfacing with existing C routines.
In order to demonstrate an application of the \verb"Matrix" class,
it was used to implement the backward error propagation algorithm
for multi-layer, feed-forward artificial neural networks.
\end{abstract}

\newcommand{\diag}{\mbox{diag}}

\section{Introduction}
The Array and Matrix classes permit two dimensional arrays to be treated as
objects that can be included in arithmetic expressions with infix operators
in much the same way as are scalar values.  The array elements are stored
in memory in row major order.  Row (column) vectors are two dimensional
arrays with just one row (column).  A matrix with exactly one row and one
column is a scalar.

\subsection{class Array}
An Array object contains four data members:
\begin{enumerate}
  \item \verb"int L" is the length of rows,
  \item \verb"int M" is the number of rows,
  \item \verb"int N" is the number of columns and
  \item \verb"double* X" is a pointer to the array.
\end{enumerate}
An array object points to a contiguous block of memory
organized into $M$ rows of $L$ numbers\footnote{The Array class may actually
be derived from single or double precision real or complex numbers.}
which represents an $M\!\times\!N$ submatrix with $N \leq L$.
See figure~\ref{fig:array}.
\begin{figure}
\newcommand{\twentyninerm}{\Large}
\input{submatrix.tex}\caption{An array object points to a contiguous block
of memory organized into $M$ rows of $L$ numbers which represents
an $M\!\times\!N$ submatrix with $N \leq L$.}\label{fig:array}
\end{figure}
Whenever an array object is declared explicitly,
the user is responsible for allocating and deallocating memory for it.
Almost all matrix operations are defined on class Array.
But many operations return an object of class Matrix.

\subsection{class Matrix}
The derived class Matrix inherits all of the data members and operations
of the base class Array.  But it automatically allocates and deallocates
memory for the array.  Storage allocated for local variables
(or intermediate results during expression evaluation) is deallocated
as soon as they exit scope.

\subsection{Indexing}
Array and Matrix class objects are indexed in exactly the same way
as two dimensional C arrays. If \verb"A" is an $M\!\times\!N$ array,
\verb"A[i]" returns a pointer of type \verb"double*"
to the first element of row~$i$ and \verb"A[i][j]" returns a reference
of type \verb"double&" to the element in column~$j$ of row~$i$
where $0\leq i<M$ and $0\leq j<N$.

\subsection{Construction}
The declaration \verb"Matrix A(m, n)" creates an uninitialized
$m\!\times\!n$ matrix.
An optional third argument, scalar $s$, in
the declaration \verb"Matrix A(m, n, s)" initializes all the elements to $s$.
The declaration \verb"Matrix x(n)" is equivalent to
the declaration \verb"Matrix x(1, n)" and creates a uninitialized
$1\!\times\!n$ row vector.
The declaration \verb"Matrix B(A)" creates a new matrix $B$ which is the
same size as matrix $A$ but allocates storage at a different location and
copies the array.  An optional third argument, pointer $p$, in
the declaration \verb"Matrix A(m, n, p)"  creates an object which points
to an $m\!\times\!n$ array at location $p$.  This declaration should be
used with care as the compiler will automatically generate code to delete
the array at location $p$ as soon as matrix $A$ exits scope.  The declaration
\verb"Array A(m, n, p)" should be used instead if the array at location $p$
must survive matrix $A$.

\subsection{Arithmetic Operations}
Unary plus, \verb"+A" $= A$, unary minus, \verb"-A" $= -A$,
scalar multiplication, \verb"s*A" $= sA = As =$ \verb"A*s",
scalar division, \verb"A/s" $= A/s$, and
element by element addition, \verb"A + B" $= A + B$,
subtraction, \verb"A - B" $= A - B$,
multiplication, \verb"A*B" $= A\!\cdot\!B$, and
division, \verb"A/B" $= A/B$, all have the expected meanings.
But other useful operations have been defined as well.

\subsubsection{Matrix Multiplication}
Matrix multiplication is effected by the infix operator, \verb"%", which
computes an inner product, \verb"A%B" $= AB^{t}$, on the rows of the operands.
This operation has indexing and memory referencing advantages\footnote{
Striding across rows when multiplying large matrices results in frequent
cache misses and page faults.} that tend to outweigh the cost of transposing
either or both of the operands.

\subsubsection{Scalar Operations}
The same scalar operation is applied to every element of the array.
\paragraph{Scalar by Matrix}
\subparagraph{Addition} \verb"s + A" $= s\vec{1}^{t}\vec{1} + A$.
\subparagraph{Subtraction} \verb"s - A" $= s\vec{1}^{t}\vec{1} - A$.
\subparagraph{Division} \verb"s/A" $= s\vec{1}^{t}\vec{1}/A$.

\paragraph{Matrix by Scalar}
\subparagraph{Addition} \verb"A + s" $= A + s\vec{1}^{t}\vec{1}$.
\subparagraph{Subtraction} \verb"A - s" $= A - s\vec{1}^{t}\vec{1}$.

\subsubsection{Vector Operations}
Element by element operations between operands of different sizes are
fatal errors unless one of the operands is a row (column) vector with
the same number of columns (rows) as the other operand\footnote{An
operand with exactly one row and one column is treated as a scalar}.
The row (column) vector is simply combined with each of
the rows (columns) of the matrix.

\paragraph{Row Vectors}
\subparagraph{Addition}
\verb"x + A" $= \vec{1}^{t}\vec{x} + A = A + \vec{1}^{t}\vec{x} =$ \verb"A + x".
\subparagraph{Subtraction}
\verb"x - A" $= \vec{1}^{t}\vec{x} - A$ and
\verb"A - x" $= A - \vec{1}^{t}\vec{x}$.
\subparagraph{Multiplication}
\verb"x*A" $= A\diag(\vec{x}) =$ \verb"A*x".
\subparagraph{Division}
\verb"x/A" $= \vec{1}^{t}\vec{x}/A$ and \verb"A/x" $= A/(\vec{1}^{t}\vec{x})$.

\paragraph{Column Vectors}
\subparagraph{Addition}
\verb"x + A" $= \vec{x}\vec{1} + A = A + \vec{x}\vec{1} =$ \verb"A + x".
\subparagraph{Subtraction}
\verb"x - A" $= \vec{x}\vec{1} - A$ and
\verb"A - x" $= A - \vec{x}\vec{1}$.
\subparagraph{Multiplication}
\verb"x*A" $= \diag(\vec{x})A =$ \verb"A*x".
\subparagraph{Division}
\verb"x/A" $= \vec{x}\vec{1}/A$ and \verb"A/x" $= A/(\vec{x}\vec{1})$.

\subsection{Shift Operators}
If the left operand is an array, $A$, and the right operand is an integer, $n$,
the shift operator \verb"<<" (\verb">>") returns a copy of $A$
shifted $n$ columns left (right).
The operator \verb"<<" fills the empty columns with zeros.
The operator \verb">>" copies the leftmost column into the empty columns.
If the left operand is an \verb"istream" and the right operand is an array,
the right shift operator, \verb">>", inputs numbers separated by white spaces
in row major order until the entire array is filled.
If the left operand is an \verb"ostream" and the right operand is an array,
the left shift operator, \verb"<<", outputs the entire array in row major order
using the default output format and terminating each row with a new line.
The \verb"format" function specifies the default output format
and returns an empty character string.  It has three arguments.
The first argument specifies the output format (initially \verb'"%g"')
for each element of the array.
The second argument (optional) specifies the number (initially \verb'4')
of elements to display on each line.
The third argument (optional) specifies the character string
(initially \verb'" "') which separate each element displayed on a line.

\subsection{Other Operators}
Two arrays with the same number of columns can be stacked,
\verb"(A^B)" $= \left[\begin{array}{c}A\\B\end{array}\right]$.
Two arrays with the same number of rows can be adjoined,
\verb"(A|B)" $= \left[A|B\right]$.
The outer product, \verb"(A&B)" $= A^{t}B$, on the rows of the operands
should be avoided unless both operands are row vectors.

\subsection{Assignment Operators}
The assignment operator, \verb"=", copies each element of its right hand side
into the corresponding element on its left hand side if both matrices are the
same size.  If the right hand side is a row (column) vector with the same
number of rows (columns) of the matrix on the left hand side, then it is copied
into each of the rows (columns) of the matrix.  If the right hand side is a
scalar, then it is copied into every element of the matrix on the left hand
side.  Operators \verb"+=", \verb"-=", \verb"*=", \verb"/=", \verb"<<="
and \verb">>=" were implemented with the expected meaning.
But operator \verb"%=" was not implemented.

\subsection{Member Functions}
\subsubsection{Data Members}
Although access to the data members, \verb"L", \verb"M", \verb"N" and \verb"X",
is not restricted in the implementation, the user must admit the possibility
of an alternate implementation and should reference them only through the
respective member functions, \verb"l()", \verb"m()", \verb"n()" and \verb"x()".

\subsubsection{Sub Array}
\verb"A.s(i, m, j, n)" is an $m\!\times\!n$ sub array
beginning with row~$i$ and column~$j$ of array (matrix) $A$.
\verb"A.s(i, m, j)" is an $m\!\times\!(N\!-\!j)$ sub array
beginning with row~$i$ and column~$j$.
\verb"A.s(i, m)" is $m$ rows beginning with row~$i$.
\verb"A.s(i)" is row~$i$.  And \verb"A.s()" is just row~$0$.

\subsubsection{Sub Matrix}
Since the sub array member function is an array object,
it should only be applied to Array or Matrix objects
declared explicitly by the user.  It should never
be applied to intermediate results in expressions
because the storage allocated for the intermediate result
may be reallocated upon return from the sub array function.
The sub matrix member function is similar to the sub array function
except that it actually allocates storage and copies the sub array into it.
The intermediate result can safely be destroyed since it is no longer needed.

\verb"A._(i, m, j, n)" returns an $m\!\times\!n$ sub matrix
beginning with row~$i$ and column~$j$ of array (matrix) $A$.
\verb"A._(i, m, j)" returns an $m\!\times\!(N\!-\!j)$ sub matrix
beginning with row~$i$ and column~$j$.
\verb"A._(i, m)" returns $m$ rows beginning with row~$i$.
\verb"A._(i)" is row~$i$.  And \verb"A._()" returns just row~$0$.

\subsubsection{Transpose}
\verb"A.t()" $= A^{t}$.

\subsubsection{Inverse}
\verb"A.i()" $= A^{-1}$.

\subsubsection{Row Sum}
\verb"A.sum()" $= \vec{1}A^{t}$.

\subsubsection{Row Sum of Squares}
\verb"A.sumsq()" $= \vec{1}(A\!\cdot\!A)^{t}$.

\subsubsection{Mapping}
\verb"A.map(f)" $= [f(a_{ij})]$.

\subsubsection{Row Minimums}
\verb"A.min()" $= [\min_{j}\{a_{ij}\}]$.

\subsubsection{Row Maximums}
\verb"A.max()" $= [\max_{j}\{a_{ij}\}]$.

\subsubsection{Minimum Index}
\verb"A.min_index()" $= iL\!+\!j$ where $a_{ij} = \min_{i,j}\{a_{ij}\}$.

\subsubsection{Maximum Index}
\verb"A.max_index()" $= iL\!+\!j$ where $a_{ij} = \max_{i,j}\{a_{ij}\}$.

\subsubsection{Index}
\verb"A.index(s)" $= iL\!+\!j$ where $a_{ij} = s$.

\subsubsection{Resize}
\verb"A.resize(m, n)" converts $A$ to an $m\!\times\!n$ matrix.
But \verb"resize" can be called with any of the arguments passed
to Matrix constructors.

\subsection{Comparison Operators}
A comparison is always valid when one of the operands in a comparison
is a scalar.  For example, the expression \verb"s == A" is true
if and only if $s=a_{ij}$ for all $i$ and $j$.
A comparison of the form \verb"A @ B"
where \verb"@" $\in \{$\verb"<", \verb"<=", \verb"==", \verb">=", \verb">"$\}$
is valid if and only if \verb"B - A" is valid
and true if and only if \verb"0 @ B - A" is true.

\section{Basetypes}
The Matrix class was implemented for \verb"double" and \verb"complex" basetypes.
The fast fourier cosine transform, \verb"ffct", and the signum function,
\verb"sgn", are included in the \verb"doubleArray" class.
The \verb"complexArray" class has no \verb"<", \verb"<=", \verb">"
or \verb">=" operators and no \verb"min", \verb"max", \verb"min_index"
or \verb"max_index" member functions but includes \verb"fft", \verb"polar",
\verb"conj", \verb"real", \verb"imag", \verb"norm" and \verb"arg" functions.
Both classes include \verb"exp", \verb"log", \verb"sin", \verb"cos", \verb"tan",
\verb"sinh", \verb"cosh", \verb"tanh", \verb"sqrt" and \verb"abs" functions.

\section{Genclass}
Instead of maintaining separate sources for both basetypes,
the Gnu container class prototype mechanism\cite{lea} was used.
The \verb"Matrix.hP" header prototype file contains the class definitions
and short inline functions.
The \verb"Matrix.ccP" source prototype file contains code
for the longer library routines.
The Gnu \verb"genclass" utility replaces every occurrence of \verb"<T>"
in these files with the desired basetype.  But since the basetype is not
the only difference between the container classes, the \verb"genclass"
shell script was modified to use the \verb"m4" macro processor
which selectively includes code depending upon the basetype.  The\\
\verb"genclass double  val Matrix" or\\
\verb"genclass complex val Matrix"
command generates the respective header and source files,
\verb"double.Matrix.h" and \verb"double.Matrix.C" or
\verb"complex.Matrix.h" and \verb"complex.Matrix.C",
for the \verb"double" or \verb"complex" basetypes respectively.

\section{Matrix Library}
The Matrix class is intended to serve as a convenient interface to existing
software.  Matrix inversion and fast fourier transforms have been implemented
with routines adapted from ``Numerical Recipes in C''\cite{press}.
The singular value decomposition algorithm, \verb"svdcmp", was used
to implement matrix inversion.  The \verb"four1", \verb"realft" and
\verb"cosft" routines were used to implement fourier transforms.
The \verb"doubleMatrix.C", \verb"complexMatrix.C", \verb"four1.C",
\verb"realft.C" and \verb"cosft.C" source files are compiled to object
files using the \verb"g++" compiler.  Then the \verb"ar" program is used
to combine the object files into a single library archive, \verb"libMatrix.a".
And the \verb"ranlib" program is used to add a table of contents so that
the routines can be linked more rapidly.

\section{Backprop}
The backward error propagation algorithm is illustrated in figure~\ref{fig:bp}.
An input/output pair is selected at random from the training sample.
The input, $\vec{x}$, is fed into the bottom layer, $0$, of the network.
The signals, $\vec{x}_{k}$, feed forward, up through the network.
The output of each layer is the input for the next layer.
The difference between the desired output, $\vec{y}$,
and the output, $\vec{x}_{K}$, of the top layer, $K\!-\!1$,
is used to compute an equivalent error, $\vec{\delta}_{K\!-\!1}$,
which is fed back into the top layer.
The equivalent errors, $\vec{\delta}_{k}$,
propagate backward, down through the network.
The equivalent errors are used to update the biases, $\vec{b}_{k}$,
and connection weights, $W_{k}$, in each layer.
\begin{figure}
\newcommand{\twentyninerm}{\Huge}
\input{backprop.tex}%
\caption{The signals, $\vec{x}_{k}$, feed forward, up through the layers.
The equivalent errors, $\vec{\delta}_{k}$, feed backward,
down through the layers.
The equivalent errors are used to update the biases, $\vec{b}_{k}$,
and the connection weights, $W_{k}$.}\label{fig:bp}
\end{figure}

The backward error propagation (backprop) algorithm
can be summarized in three steps:
\begin{enumerate}
\item Feed Forward
\begin{equation}
\vec{x}_{0} = \vec{x}
\end{equation}
\begin{equation}
\vec{x}_{k\!+\!1} = f(\vec{b}_{k} + \vec{x}_{k}W_{k}^{t})
\end{equation}
\begin{equation}
\vec{x}_{k\!+\!1}^{\prime} = f^{\prime}(\vec{b}_{k} + \vec{x}_{k}W_{k}^{t})
\end{equation}
\item Back Propagate
\begin{equation}
\vec{\delta}_{K\!-\!1} = \eta(\vec{y} - \vec{x}_{K})\cdot\vec{x}_{K}^{\prime}
\end{equation}
\begin{equation}
\vec{\delta}_{k\!-\!1} = \vec{\delta}_{k}W_{k}\cdot\vec{x}_{k}^{\prime}
\end{equation}
\item Update
\begin{equation}
\vec{b}_{k} \leftarrow \vec{b}_{k} + \vec{\delta}_{k}
\end{equation}
\begin{equation}
W_{k} \leftarrow W_{k} + \vec{\delta}_{k}^{t}\vec{x}_{k}
\end{equation}
\end{enumerate}
where the learning rate, $\eta$, is a small number.

A multi-layer, feed-forward, artificial neural network
can approximate any continuous function arbitrarily closely
provided that there are enough (non-linear) neurons in the hidden layers.
The output neurons are linear and the hyperbolic tangent, $\tanh(\cdot)$,
is used for the non-linear hidden neurons so that the derivatives,
$\vec{x}_{k}^\prime = (\vec{1} - \vec{x}_{k}\cdot\vec{x}_{k})$,
are easily calculated.

The backprop algorithm itself is easy to implement with the C++ Matrix class.
\begin{verbatim}
x[0] = x;
for (k = 0; k < K-1; k++)			// Feed Forward
  x[k+1] = tanh(b[k] + x[k]%W[k]);
x[K] = b[K-1] + x[K-1]%W[K-1];

d[K-1] = eta*(y - x[K]);			// Back Propagate
for (k = K-1; k > 0; k--)
  d[k-1] = d[k]%W[k].t()*(1.0 - x[k]*x[k]);

for (k = 0; k < K; k++)				// Update
  {
    b[k] += d[k];
    W[k] += d[k]&x[k];
  };
\end{verbatim}
The rest of the \verb"backprop" program processes command line options,
constructs arrays of matrices for the network and the training sample
and reads data into them then writes the network matrices out again
before it terminates.
Network and training sample data are kept in separate files
with different formats detailed in table~\ref{tab:nff}
and table~\ref{tab:sff} respectively.
\begin{table}
\begin{tabular}{ll}
$K$					& \# of layers\\
$n_{0}$					& \# of inputs to layer $0$\\
$n_{1}$					& \# of inputs to layer $1$\\
$[\vec{b}_{0}^{t}, W_{0}]$		& layer $0$ biases \& weights\\
$n_{2}$					& \# of inputs to layer $2$\\
$[\vec{b}_{1}^{t}, W_{1}]$		& layer $1$ biases \& weights\\
$\vdots$\\
$n_{K}$					& number of outputs\\
$[\vec{b}_{K\!-\!1}^{t}, W_{K\!-\!1}]$	& layer $K\!-\!1$ biases \& weights\\
\end{tabular}\caption{Network file format.}\label{tab:nff}
\end{table}
\begin{table}
\begin{tabular}{ll}
$m$					& \# of examples\\
$(\vec{x}_{0}, \vec{y}_{0})$		& example $0$\\
$(\vec{x}_{1}, \vec{y}_{1})$		& example $1$\\
$\vdots$\\
$(\vec{x}_{m\!-\!1}, \vec{y}_{m\!-\!1})$& example $m\!-\!1$\\
\end{tabular}\caption{Sample file format.}\label{tab:sff}
\end{table}

\subsection{Spiral Problem}
A three layer network, $N(2, 5, 10, 2)$, with two inputs,
five hidden neurons in the first layer, ten hidden neurons in the second layer
and two output neurons was used trained to compute the Cartesian coordinates,
\begin{equation}
(x, y)	= \vec{f}(\rho, \phi)
	= (\rho\cos(\rho - \phi), \rho\sin(\rho - \phi)),
\end{equation}
of a point on a spiral where $0 \leq \rho < \pi$ and  $0 \leq \phi < \pi$
are the radius and angle respectively.

The network biases and connection weights were initialized
to small random numbers and stored in file \verb"ffnet.new".
The training sample is stored in file \verb"ffnet.S".
Training is accomplished in three steps:
\begin{enumerate}
\item \verb"mv ffnet.new ffnet.old"\\
      Move the file \verb"ffnet.new" to file \verb"ffnet.old".
\item \verb"cat ffnet.old ffnet.S | (backprop -v > ffnet.new) >>& ffnet.err &"\\
      The program reads from standard input and writes to standard output.
      The \verb"ffnet.old" and \verb"ffnet.S" files are concatenate together
      and piped to the \verb"backprop" program.  The \verb"backprop" output 
      is redirected to file \verb"ffnet.new".  And any error messages are
      appended to file \verb"ffnet.err".  It runs in the background so that
      the user can continue working until the program terminates.
\item \verb"cat ffnet.new ffnet.S | evaluate | graph -m 0 | plot"\\
      In order to evaluate the network performance,
      the user plots the network output versus the desired output
      for every example in the training sample.
      The \verb"ffnet.new" and \verb"ffnet.S" files are concatenated together,
      piped through \verb"evaluate" and \verb"graph" to \verb"plot" which
      displays the plot on the user's terminal.
\end{enumerate}
The previous steps can be repeated until the user concludes
that the performance is acceptable
or that the network will not learn the function.

\appendix
\section{Notation}
An $N$ element row or column vector is usually denoted by
a lower case symbol (i.e. $\vec{a}$). And an $M\!\times\!N$
matrix is usually denoted by an upper case symbol (i.e. $A$).
Row $0\leq i<M$ of matrix $A$ is denoted $\vec{a}_{i}$.
Column $0\leq j<N$ of matrix $A$ is denoted $\vec{a}_{j}$.
And $a_{ij}$ denotes the element in column~$j$ of row~$i$.
The inverse and transpose of matrix $A$ are denoted $A^{-1}$ and
$A^{t}$ respectively.  The vector $\vec{1}$ is a row vector of $1$s.
The matrix $\diag(\vec{a})$ is an $N\!\times\!N$ diagonal matrix
where the diagonal elements are the $N$ elements of vector $\vec{a}$.
The expression $[a_{j}]$ is a row vector and $[a_{ij}]$ is a matrix
composed of scalar elements, $a_{j}$ and $a_{ij}$ respectively.
But $[\vec{a}_{j}]$ is a matrix composed of column vectors, $\vec{a}_{j}$.
Matrix multiplication, $AB$, and scalar multiplication, $sA = As$,
are denoted by juxtaposition.  But scalar division is denoted, $A/s$.
Element by element addition, subtraction, multiplication, division and
function evaluation are denoted $A + B$, $A - B$, $A\!\cdot\!B$, $A/B$
and $f(A)$ respectively.

\nocite{stroustrup}
\bibliography{Matrix}
\end{document}
